%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FUNCTION IIV_cereal_code_example_usage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DESCRIPTION: This function applies the IIV_Bounds function to the
%   Cereal dataset used in Nevo and Rosen (2012) [NR]. The files 
%   Cereal_Restat_citycode.csv and IIV_Bounds.m must be in
%   the same directory.
% INPUT PARAMETERS:
%   R: number of simulation draws for computing critical values.
%   bootstrap_samples: number of draws to use when bootstrapping.
% OUTPUT PARAMETERS:
%   Bound estimates and confidence intervals for slope coefficients on price
%   and advertising under restrictions A3 (col4) and restrictions
%   A3 and A4 (col5) for the application described in NR Section V.
% EXAMPLE USAGE:
%   At the Matlab command prompt execute:
%   [col4_bound_estimates,col5_bound_estimates,col4_price_CI,...
%   col4_adv_CI,col5_price_CI,col5_adv_CI] = ...
%   IIV_cereal_code_example_usage(1000,10000)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [col4_bound_estimates,col5_bound_estimates,col4_price_CI,...
    col4_adv_CI,col5_price_CI,col5_adv_CI] = ...
    IIV_cereal_code_example_usage(bootstrap_samples,R)

% LOAD DATA
data = csvread('Cereal_Restat_citycode.csv',1,0);
n = size(data,1);

% SET VARIABLES
y = data(:,52);
price = data(:,3);
data(:,6) = data(:,6)/10; % put in adv/10 as in Table 2
city = data(:,58); % 1 if SF, 0 if BOS
qavgp = data(:,53);
p_SF = data(:,54);
p_BS = data(:,55);
qavgpo = city .* p_BS + (1-city) .* p_SF;
t1 = std(qavgp);
t2 = std(qavgpo);
ivdif = qavgpo - qavgp;
ivdif1 = (t1*qavgpo - t2*qavgp)/(t1+t2); % Not used, corresponds to alternative gamma, as discussed in paper

% For use in computing additional upper bounds in Proposition 5
sz = std(ivdif);
sx = std(price);
divdif = sx*ivdif - sz*price; % divdif is vstar1 in Proposition 5
v1_1 = sx*qavgp - std(qavgp)*price; % V1(lambda) with lambda = 1
v2_1 = sx*qavgpo - std(qavgpo)*price; % V2(lambda) with lambda = 1
 
regvars = data(:,[6:29,31,33:51,58]); % bd24 ommitted
X_OLS = [price,regvars,ones(n,1)];
Z1 = [qavgp,regvars,ones(n,1)];
Z2 = [qavgpo,regvars,ones(n,1)];
Zivdif = [ivdif,regvars,ones(n,1)];
Zdivdif = [divdif,regvars,ones(n,1)];
Zv1_1 = [v1_1,regvars,ones(n,1)];
Zv2_1 = [v2_1,regvars,ones(n,1)];

%------------ RUN THE REGRESSIONS TO PRODUCE TABLE 2 BOUND ESTIMATES

[beta_OLS,var_OLS] = TWOSLS_REG(y,X_OLS,X_OLS); % Upper bound from first line of Prop 5
se_OLS = sqrt(diag(var_OLS));

[beta_qavgp,var_qavgp] = TWOSLS_REG(y,X_OLS,Z1); % Upper bound from both lines of Prop 5
se_qavgp = sqrt(diag(var_qavgp));

[beta_qavgpo,var_qavgpo] = TWOSLS_REG(y,X_OLS,Z2); % Upper bound from both lines of Prop 5
se_qavgpo = sqrt(diag(var_qavgpo));

[beta_ivdif,var_ivdif] = TWOSLS_REG(y,X_OLS,Zivdif); % This is the lower bound for both column 4 and column 5 (both parts of Prop 5)
se_ivdif = sqrt(diag(var_ivdif));

[beta_divdif,var_divdif] = TWOSLS_REG(y,X_OLS,Zdivdif); % Upper bound from second line of Prop 5
se_divdif = sqrt(diag(var_divdif));

[beta_v1_1,var_v1_1] = TWOSLS_REG(y,X_OLS,Zv1_1); % Upper bound from second line of Prop 5
se_v1_1 = sqrt(diag(var_v1_1));

[beta_v2_1,var_v2_1] = TWOSLS_REG(y,X_OLS,Zv2_1); % Upper bound from second line of Prop 5
se_v2_1 = sqrt(diag(var_v2_1));

%----------------------COLUMN 4-------------------------------------------- 
% Collect column 4 bound estimates using first results in Proposition 5

upper_bound_price_estimates_column4 = [beta_OLS(1),beta_qavgp(1),beta_qavgpo(1)]';
lower_bound_price_estimates_column4 = beta_ivdif(1);
upper_bound_price_se_column4 = [se_OLS(1),se_qavgp(1),se_qavgpo(1)];
lower_bound_price_se_column4 = se_ivdif(1);

upper_bound_adv_estimates_column4 = [beta_OLS(2),beta_qavgp(2),beta_qavgpo(2)]';
lower_bound_adv_estimates_column4 = beta_ivdif(2);
upper_bound_adv_se_column4 = [se_OLS(2),se_qavgp(2),se_qavgpo(2)];
lower_bound_adv_se_column4 = se_ivdif(2);

%----------------------COLUMN 5-------------------------------------------- 
% Collect column 5 bound estimates using second results in Proposition 5 using A4
upper_bound_price_estimates_column5 = [beta_qavgp(1),beta_qavgpo(1),beta_v1_1(1),beta_v2_1(1),beta_divdif(1)]';
lower_bound_price_estimates_column5 = beta_ivdif(1);
upper_bound_price_se_column5 = [se_qavgp(1),se_qavgpo(1),se_v1_1(1),se_v2_1(1),se_divdif(1)];
lower_bound_price_se_column5 = se_ivdif(1);

upper_bound_adv_estimates_column5 = [beta_qavgp(2),beta_qavgpo(2),beta_v1_1(2),beta_v2_1(2),beta_divdif(2)]';
lower_bound_adv_estimates_column5 = beta_ivdif(2);
upper_bound_adv_se_column5 = [se_qavgp(2),se_qavgpo(2),se_v1_1(2),se_v2_1(2),se_divdif(2)];
lower_bound_adv_se_column5 = se_ivdif(2);

% USE BOOSTRAP RESAMPLING TO ESTIMATE THE CORRELATION MATRIX FOR THE BOUND ESTIMATES
bootstrap_draws = unidrnd(n,n,bootstrap_samples);
bootstrap_pricecoeff_estimates = zeros(bootstrap_samples,6); % Upper bounds on price coefficients
bootstrap_advcoeff_estimates = zeros(bootstrap_samples,6); % Upper bounds on adv/10 coefficients

for b=1:bootstrap_samples
   % Collect bootstrap sample b draws
   y_b = y(bootstrap_draws(:,b));
   X_OLS_b = X_OLS(bootstrap_draws(:,b),:);
   qavgp_b = qavgp(bootstrap_draws(:,b));
   qavgpo_b = qavgpo(bootstrap_draws(:,b));
     
   % Form regression variables including constant
   regvars_cons_b = X_OLS_b(:,2:size(X_OLS,2));

   % Compute variables and IVs
   ivdif_b = qavgpo_b - qavgp_b;
   price_b = X_OLS_b(:,1);
   sz_b = std(ivdif_b);
   sx_b = std(price_b);
   
   divdif_b = sx_b*ivdif_b - sz_b*price_b; % divdif is vstar1 in Proposition 5
   v1_1_b = sx_b*qavgp_b - std(qavgp_b)*price_b; % V1(lambda) with lambda = 1
   v2_1_b = sx_b*qavgpo_b - std(qavgpo_b)*price_b; % V2(lambda) with lambda = 1
   
   Z1_b = [qavgp_b,regvars_cons_b];
   Z2_b = [qavgpo_b,regvars_cons_b];
   Zdivdif_b = [divdif_b,regvars_cons_b];
   Zv1_1_b = [v1_1_b,regvars_cons_b];
   Zv2_1_b = [v2_1_b,regvars_cons_b];

   beta_price_OLS_b = (X_OLS_b'*X_OLS_b)\(X_OLS_b'*y_b);
   beta_qavgp_b = (Z1_b' * X_OLS_b) \ (Z1_b' * y_b);
   beta_qavgpo_b = (Z2_b' * X_OLS_b) \ (Z2_b' * y_b);
   beta_v1_1_b = (Zv1_1_b' * X_OLS_b) \ (Zv1_1_b' * y_b);
   beta_v2_1_b = (Zv2_1_b' * X_OLS_b) \ (Zv2_1_b' * y_b);
   beta_divdif_b = (Zdivdif_b' * X_OLS_b) \ (Zdivdif_b' * y_b);

   bootstrap_pricecoeff_estimates(b,:) = [beta_price_OLS_b(1),beta_qavgp_b(1),beta_qavgpo_b(1),beta_v1_1_b(1), beta_v2_1_b(1),beta_divdif_b(1)];
   bootstrap_advcoeff_estimates(b,:) = [beta_price_OLS_b(2),beta_qavgp_b(2),beta_qavgpo_b(2),beta_v1_1_b(2), beta_v2_1_b(2),beta_divdif_b(2)];
   
end

% Compute correlation matrices for upper bound estimates
% -- Note these are CORRELATIION matrices, not VARIANCE COVARIANCE matrices
% as indicated in step 2 on page 666. -- 
Omega_price = corr(bootstrap_pricecoeff_estimates);
Omega_adv = corr(bootstrap_advcoeff_estimates);

% Use inequality_selection = 'AIS' to produce confidence intervals using
% the adaptive inequality selection  procedure from
% Chernzhukov, Lee, and Rosen (2013).
% ** This is the recommended method. **

inequality_selection = 'AIS';

% Uncomment the following line to instead select the inequalities
% that come within c * sqrt(log(n)) * max(se) of the maximal/minimal inequality,
% where c is a fixed positive constant.

% inequality_selection = 'IIV';

% c = 2 in procedure described in Section IV of NR.
% c not used if method is AIS.

c = 2;

% c = 0.01 produces the confidence intervals of Section V of NR in Table 2.
% c = 2 is recommended if IIV is used, though that will produce slightly larger
% confidence intervals.

% -- COLUMN 4 -- BOUNDS ON PRICE COEFFICIENT ----------------------%
Omega_u = Omega_price(1:3,1:3);
se_u = [se_OLS(1),se_qavgp(1),se_qavgpo(1)]'; % standard errors
upper_bound_estimates = [beta_OLS(1),beta_qavgp(1),beta_qavgpo(1)]';
[col4_price_bounds,col4_price_CI] = IIV_Bounds(upper_bound_estimates,...
                                                  lower_bound_price_estimates_column4,...
                                                  se_u,lower_bound_price_se_column4,...
                                                  Omega_u,1,0.975,n,c,R,inequality_selection);

% -- COLUMN 4 -- BOUNDS ON ADVERTISING/10 COEFFICIENT ----------------------%
Omega_u = Omega_adv(1:3,1:3);
se_u = [se_OLS(2),se_qavgp(2),se_qavgpo(2)]'; % standard errors
upper_bound_estimates = [beta_OLS(2),beta_qavgp(2),beta_qavgpo(2)]';
[col4_adv_bounds,col4_adv_CI] = IIV_Bounds(upper_bound_estimates,...
                                                  lower_bound_adv_estimates_column4,...
                                                  se_u,lower_bound_adv_se_column4,...
                                                  Omega_u,1,0.975,n,c,R,inequality_selection);

col4_bound_estimates = [col4_price_bounds;col4_adv_bounds];

% -- COLUMN 5 -- BOUNDS ON PRICE COEFFICIENT ----------------------%
Omega_u = Omega_price(2:6,2:6);
se_u = [se_qavgp(1),se_qavgpo(1),se_divdif(1),se_v1_1(1),se_v2_1(1)]'; % standard errors
upper_bound_estimates = [beta_qavgp(1),beta_qavgpo(1),beta_divdif(1),beta_v1_1(1),beta_v2_1(1)]';
[col5_price_bounds,col5_price_CI] = IIV_Bounds(upper_bound_estimates,...
                                                  lower_bound_price_estimates_column5,...
                                                  se_u,lower_bound_price_se_column5,...
                                                  Omega_u,1,0.975,n,c,R,inequality_selection);
            
% -- COLUMN 5 -- BOUNDS ON ADVERTISING/10 COEFFICIENT ----------------------%
Omega_u = Omega_adv(2:6,2:6);
se_u = [se_qavgp(2),se_qavgpo(2),se_divdif(2),se_v1_1(2),se_v2_1(2)]'; % standard errors
upper_bound_estimates = [beta_qavgp(2),beta_qavgpo(2),beta_divdif(2),beta_v1_1(2),beta_v2_1(2)]';

[col5_adv_bounds,col5_adv_CI] = IIV_Bounds(upper_bound_estimates,...
                                                  lower_bound_adv_estimates_column5,...
                                                  se_u,lower_bound_adv_se_column5,...
                                                  Omega_u,1,0.975,n,c,R,inequality_selection);

col5_bound_estimates = [col5_price_bounds;col5_adv_bounds];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function TWOSLS_REG implements two stage least squares
% y univariate outcome variable
% X multivariate "RHS" variables
% W instruments for X
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [betahat,varcovhat] = TWOSLS_REG(y, X, W)

T = size(y,1);
k = size(X,2);

Pw=W*((W'*W)\W'); 
betahat = (X'*Pw*X)\(X'*Pw*y);

dist = (y-X*betahat)'*(y-X*betahat);
varcovhat = (dist/(T-k))*inv(X'*Pw*X);

end
